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Abstract 

Several recently developed multisymplectic schemes for Hamiltonian PDEs have been shown 
to preserve associated local conservation laws and constraints very well in long time numerical 
simulations. Backward error analysis for PDEs, or the method of modified equations, is a useful 
technique for studying the qualitative behavior of a discretization and provides insight into the 
preservation properties of the scheme. In this paper we initiate a backward error analysis for 
PDE discretizations, in particular of multisymplectic box schemes for the nonlinear Schrodinger 
equation. We show that the associated modified differential equations are also multisymplectic 
and derive the modified conservation laws which are satisfied to higher order by the numerical 
solution. Higher order preservation of the modified local conservation laws is verified numerically. 

1 Introduction 

When developing numerical integrators for Hamiltonian PDEs that possess a multisymplectic struc- 
ture (i.e. symplectic in both space and time), it is natural to require the numerical scheme to 
preserve exactly a discrete version of the multisymplectic conservation law (MSCL) [21 El- How- 
ever, this does not imply preservation of other dynamical invariants of the system such as the local 
energy and momentum conservation laws or global invariants which determine the phase space 
structure. A question that immediately arises then is, to what extent are the other invariants of 
the system preserved? Recent numerical experiments using multisymplectic integrators for non- 
linear wave equations (e.g. the nonlinear Schrodinger (NLS), sine-Gordon, and Gross-Pitaevskii 
eqautions) show that the local conservations laws are preserved very well, although not exactly, 
over long times [SJinilZI- Further, the improved preservation of the local conservation laws is re- 
flected in an improved preservation of complicated phase space structures 7 . This is reminiscent 
of the behavior of symplectic schemes for Hamiltonian ODEs. Symplectic integrators are designed 
to preserve the symplectic structure, not to preserve the energy. In fact, for general Hamiltonian 
systems, conservation of the symplectic structure and conservation of energy are conflicting require- 
ments that, in general, are not solved simultaneously by a given scheme @J. Even so, symplectic 
integrators preserve the Hamiltonian extremely well over very long times. 

Backward error analysis (BEA), or the method of modified equations, is a particularly insightful 
technique for studying the qualitative behavior of a discretization as well as an alternative method 
for checking the accuracy of the numerical solution ^2] • Since our main interest lies in the geometry 
preserving properties of multisymplectic schemes, the main question backward error analysis tries 
to answer (whether the distinguishing properties of the original equation carry over to the modified 
equation which the numerical solution satisfies to higher order) becomes relevant to our study. For 
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a given scheme, the derivation of the associated modified equation is related to the calculation of 
the local truncation error and has, typically, been used to examine the dispersive, dissipative and 
diffusive properties of PDE discretizations. For example, in the numerical analysis of linear PDEs 
a backward error analysis of the Lax-Friedrichs method or the upwind method for the advection 
equation produces in both cases a modified equation that is an advection-diffusion equation. This 
helps one to understand the qualitative behavior of the methods and, from this perspective, explains 
why the numerical solution in both cases becomes smeared out as time evolves. 

Likewise, BEA is an important tool in the study of geometric integrators @J I1UI [5J. For 
Hamiltonian ODEs, symplectic methods lead to modified equations which are also Hamiltonian. In 
fact, the modified equation of a Hamiltonian ODE is also Hamiltonian if and only if the integrator 
is symplectic; this is then used to rigorously establish that a symplectic integrator almost preserves 
the total energy over an exponentially long period of time 4 . In striking contrast, nonsymplectic 
methods used to integrate Hamiltonian ODEs can introduce dissipation, a feature which is readily 
predicted by the dissipative form of the modified equations. Less has been established using BEA 
for Hamiltonian PDEs since there are a variety of ways to implement a BEA and the relevance of 
the analysis is open to interpretation. Spatial discretization of a PDE results in a system of ODES 
to which a standard BEA can be applied to derive a modified equation that is satisfied to higher 
order in one independent variable. Alternatively, a BEA can be used to derive modified equations 
for the PDE that are satisfied to higher order in both space and time jlUl 15]. 

In this paper we implement a formal backward error analysis in both space and time of two 
multisymplectic box schemes, the Euler and the centered cell box schemes, as applied to the non- 
linear Schrodinger equation. We find that the modified equations of these box schemes are also 
multisymplectic. The modified PDEs are used to derive modified conservation laws of energy and 
momentum that are approximated by the MS scheme to higher order in space and time. For the 
centered cell discretization of the NLS we numerically verify that the modified conservation laws 
are satisfied to higher order by the numerical solution. This provides a partial explanation of the 
superior resolution of the local conservation laws and global invariants by MS schemes (e.g. see the 
numerical experiments in section 5) and a deeper understanding of the local and global properties 
of MS integrators. 

The paper is organized as follows. In the next section we recall the multisymplectic formulation 
of Hamiltonian PDEs and of the NLS equation. In section 3 we introduce the box schemes, establish 
multisymplecticity, and apply them to the NLS equation. We present a straightforward method 
for obtaining compact box schemes that is applicable to many multisymplectic PDEs. Section 4 
contains the backward error analysis of the discretizations. In section 5 numerical experiments 
for the MS centered cell box scheme are discussed, illustrating the remarkable behavior of MS 
schemes. Higher order preservation of the modified local conservation laws is verified numerically, 
which supports the use of MS integrators in long time numerical simulations of Hamiltonian PDEs. 

2 Multisymplectic Hamiltonian PDEs 

A Hamiltonian PDE (in the "1+1" case) is said to be multisymplectic if it can be written as 

Mz t + Kz x = V Z S, z G IT, (1) 

where M, K 6 M nXT1 are skew-symmetric matrices and S : JR n — * H is a smooth function of the 
state variable z(x, t) [12 Ej- The variational equation associated with Q is given by 

Mdz t + K dz x = S zz dz. (2) 
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The Hamiltonian system (^Q) is multisymplectic in the sense that associated with M and K are the 
2-forms 

uj = -(dz AMdz), K = -(dzAKdz), (3) 

which define a symplectic space-time structure (symplectic with respect to more than one indepen- 
dent variable). 

Any system of the form Q satisfies conservation of symplecticity. Let dz be any solution of 
the variational equation @. Then it can be shown that uj and k, as defined in (J3J), satisfy the 
multisymplectic conservation law (MSCL): 



duj dn 
dt dx 



0. 



(4) 



This result is obtained by noting that 
2uj t = {dz A Mdz) t = 



dz t A Mdz + dz A Mdz t 

- (Mdz t ) A dz + dz A Mdz t 

— (S zz dz — K.dz x ) A dz + dz A (S zz dz 

— (dz x A Kdz + dz A Kdz x ) 

- (dz A Kdz) x = -2k x 



Kdz x ) 



since M, K are skew-symmetric and S zz is symmetric. The MSCL (J3J is a local property and 
expresses the fact that symplecticity for Hamiltonian PDEs can vary locally over the spatial domain. 

An important consequence of the MS structure is that when the Hamiltonian S(z) is independent 
of t and x, the PDE has local energy and momentum conservation laws [111 [2] 



E t + F x 
h + G x 



0: 
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E = S(z) + h^Kz, 
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G = S(z) + ±Zf Mz, 
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-^zjKz, 
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(5) 
(6) 



For periodic boundary conditions, the local conservation laws can be integrated in x to obtain 
global conservation of energy and momentum. 



2.1 Multisymplectic formulation of the NLS equation 

The focusing one dimensional nonlinear Schrddinger (NLS) equation, 

i 1 9 

iut + u xx + 2 kt u = 0, 



(7) 



can be written in multisymplectic form by letting u = p + iq and introducing the new variables 
v = p x , w = q x . Separating Q into real and imaginary parts, we obtain the system [3]: 

v x = 2 (p 2 + q 2 ) p 

„2 , „2\ 
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which is equivalent to the multisymplectic form Q for the NLS equation with 
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and Hamiltonian 



s= l - 

2 
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Implementing for the NLS equation yields the local energy conservation law (LECL) 
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p 2 + <? 2 



2 2 
V — W 



, F = vp t +wq t , 



(9) 



3t + ^ = 0, 

and the local momentum conservation law (LMCL) 

J t + G a . = 0, I = pw-qv, G=(f + q 2 ) 2 +v 2 +w 2 -{pq t -p t q). (10) 
Additionally we have a norm conservation law for the NLS equation 



i\T t + M x = 0, AT = I (p 2 + q 2 ) , M 



qv — pw. 



(11) 



These three equations, when integrated with respect to x, yield the classic global conservation of 
energy £ (t) (Hamiltonian), momentum X{t) and norm M(t). 



3 Multisymplectic box schemes 

Multisymplectic discretizations are numerical schemes for approximating Q which preserve a dis- 
crete version of the multisymplectic conservation law (JIJ). That is, if the discretization of the 
multisymplectic PDE and its conservation law are written schematically as 



and 



M^ j 4 + Kd^4 = (y z s(4)) j . 



& <4 + au 4 = 0, 



(12) 



(13) 



where ff = f(xi,tj), and d 1 ^ and d 1 ^ are discretizations of the corresponding derivatives dt and 
d x , then the numerical scheme (|12j) is said to be multisymplectic if (|13|) is a discrete conservation 
law of {I2J) dUE]. 

A standard method for constructing multisymplectic schemes is to apply a known symplectic 
discretization to each independent variable. For example, splitting the matrices M and K as 



M = M + + ML and K = K+ + K_ with MT = -M_ and KT = -K_, 



(14) 



and using the symplectic Euler forward-backward difference approximations on both space and 
time derivatives yields the Euler box scheme 



z l z 



+ K. 



z z 



A/ A/ — A/ — A/ V * S (*%) 



(15) 



Similarly, applying the symplectic midpoint rule to both the time and space derivatives in (JJ) 
yields a "centered cell" box discretization 



At 



Ax 



1/2 
'1/2 



(16) 
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where 

Z l/2 = \ ( Z + Z l) ^ Z i = \ ( z i + ^) ' Z l/2 = 4 ( Z + Z + Z l + Z l) • (17) 

The local truncation error for the Euler box scheme is O (At + Ax 2 ) , while for the centered cell 
discretization it is O (At 2 + Ax 2 ). 

Multisymplecticity of schemes (fT5]l and (|T6|) is easily established. For example, to do so for the 
centered cell scheme, we use the discrete variational equation associated with (jlfij) given by 

(dz\,~ - dz9„\ fdzV 2 - dzl /2 \ , /, 

M (^_JZi) + K (^^-J = S. <. (IS) 

1 /2 

Taking the wedge product of dzA 2 with (|18|). note that the right-hand side is zero, since S z ^ is 
symmetric. The terms on the left-hand side can be simplified 

dz\ /2 A M (dz\ /2 - dz° 1/2 ) = \ (dz\ /2 + dz° 1/2 ) A M (dz\ /2 - dz° 1/2 ) 

= \ ( dz l/2 A Mdz{ /2 - dz° /2 A M(fe° /2 ) 



u l/2 ~ U l/2> 



whereas, 



dz\ /2 A K (dz\ 12 - dz]' 2 ) = \ (dz] 12 + dz 1 ' 2 ) A K (dz 1 ! 2 _ ^ 

= \ (dz\ /2 A Kdz{ /2 - dzl 12 A Kdzl /2 



1/2 1/2 



This implies that the numerical scheme 1)16(1 satisfies the discrete multisymplectic conservation law 

/, ,1 , ,0 \ / 1/2 1/2 \ 

3.1 Multisymplectic box schemes for the NLS equation 

The multisymplectic centered cell box scheme was first developed for the NLS equation in [2] 
where an apparently ad hoc reduction provided a particularly compact form of the scheme. This 
reduction turns out to be generalizable as can be seen in McLachlan's derivation of box schemes for 
the Korteweg de Vries equation [Q. Here we present a general approach for constructing compact 
box schemes which is applicable to many multisymplectic PDEs. 

3.1.1 Euler box scheme for the NLS equation 

We begin by introducing the following finite difference operators 

~±l o J _ J 

Dfz = ± Zi ~ Zi and £)± Z = ±%-^. 
* At x Ax 

In terms of these operators the Euler box scheme ()15() takes the form 

M + D+z + M-Dtz + K + D+z + K.D'z = V z S{z$). (19) 
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For the NLS, M and K are split using (|14|). where 
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After eliminating v and w the system reduces to 
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where we have set D 2 = D^D X = D x D x . When the second equation is shifted in time, the 
resulting six-point box scheme in stencil format is : 
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3.1.2 Centered cell box scheme for the NLS equation 

As before, we begin by introducing the appropriate finite difference operators 



Mtz 



M x z 



zi + zi 



D t z 



D x z 



2 ' 2 ' At 

In terms of these operators, the centered-cell discretization ()16|) becomes 

MD t M x z + KD x M t z = V Z S (M x M t z) , 

with discrete conservation law 

dz A MD t M x dz + dz A KD x M t dz = 0. 

The system which results upon applying (|2~Tj) to © is 

2 \{M x M t p) 2 + (M x M t q) 2 ] M x M t p 
2 [{M x M t pf + {M x M t q) 2 \ M x M t q 
M x M t v 



Ax 



(20) 
(21) 



D t M x q - D x M t v 
-D t M x p - D x M t w 

D x M t p 
D x M t q 



(22) 



M x M t w. 



Since the operators in (|20|) commute, by multiplying the first two equations in (|22[) by M x and 
back substituting v and w into the first two equations we obtain 

D t M 2 q - D 2 x M t p = 2M X ([(M x M t p) 2 + (M x M t qf] M x M t p) , 



-D t M 2 p - D 2 x M t q 



2M X ([(M x M t p) 2 + (M x M t q) 2 ] M x M t q 
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Recombining these equations into a single complex equation (with u = p + iq) yields the multisym- 
plectic box scheme for the NLS equation 



iD t M%u + D 2 x M t u - 2M X (\M x M t u\ 2 M x M t uj = 0, 



or equivalently 
,1 i „,i 
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In finite difference stencil format the six-point box scheme is given by 
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The centered cell scheme naturally gives a two time level stencil for the NLS equation. If every term 
in (J23J) contained a common factor, e.g. M x or Mt, further compactification would be possible. As 
it is, an additional reduction of (|23[) is not possible. 

4 Backward Error Analysis 

A useful method for analysing the qualitative behavior of symplectic methods for ODEs has been 
backward error analysis, where one interprets the numerical solution as the "nearly" exact so- 
lution of a modified Hamiltonian differential equation. In this section we implement a BEA in 
space and time for the multisymplectic box schemes. The modified differential equations are also 
multisymplectic and satisfy modified conservation laws. 



4.1 BEA for the Euler box scheme 

Let z be a differentiable function that, when evaluated at the lattice points, satisfies the Euler box 
scheme 1)19(1 . Using the Taylor series expansions in t about z = z(xi,tj) 

4** = z±Atz t + \At 2 zu ± • • • 
and equivalent expansions in x, we obtain to first order the following modified equation 

Mz t + \At (M+ - M_) z tt + Kz x + \Ax (K+ - K_) z xx = V z S(z). (25) 
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equation 1)25(1 can also be written in the multisymplectic form 
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and 

S(z) = S - ^AtzJ'Nzt - \Axz x Lz x . 
Applying equation (|25() to the NLS system and eliminating v and w yields the reduced system 

q t + I At q tt - p xx + \ Ax 2 p xxxx = 2 (p 2 + q 2 ) p 
-pt + \Atp u -q xx + \Ax 2 q xxxx = 2(p 2 +q 2 )q, 

or setting u = p + iq, the single equation 

iu t + u xx + 2\u\ 2 u + \Atu u + ^Ax 2 u xxxx = 0, 
which is an O (At + Ax 2 ) perturbation of the NLS. 
4.2 BEA for the centered cell box scheme 

We now assume z is a sufficiently smooth function that, when evaluated at the lattice points, is 
a solution to the centered cell scheme ([21 j). Expanding z in a Taylor series about the midpoints 
(xi+i/2,t j+1 / 2 ) we obtain 



At i f At\ 2 1 /At\ 3 

where to simplify the notation and 1 denote the grid points, 1/2 denotes the midpoints, and 
z = z(xi/2, t\/2)- The symplectic midpoint rule approximation of the time derivative is given by 

and, similarly, the space derivative is approximated by 

.1/2 .1/2 Ax 2 



Ax 24 

Substituting these expansions into (|16|). one finds that, to order 0(At 4 + Ax 4 ), z satisfies the 
modified PDE 

At 2 Ax 2 
Mz t + — M % + Kz x + —Kz xxx = V z S(z), (26) 

where all quantities are evaluated at the midpoint z = £1/2)- 

When applying equation (|2fij) to the NLS example, the resulting modified system of equations 
can be reduced to 

. l2 Ax 2 Ax 2 
itk + u xx + 2\u\ u = ~ % -^ u ttt —u xxxx . 

which is an 0(Ai 2 + Ax 2 ) perturbation of NLS. 

The modified local conservation laws can be obtained directly from equation Q26|) by multiplying 
the equation from the left by z% to obtain an energy conservation law and by z x to obtain a 
momentum conservation law. We prefer to show that the modified equation can be written in MS 
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form and from this formulation obtain the associated local conservation laws via equations 
Introducing the augmented variables 

z = (z, zt, ztt, z x ,z xx ) T , S = S + -^-z[ t Mzt + -^-z xx Kz x 

the modified equations (|26j) can be written in the MS form 

Mz t + Kz x = V~ z S(z), 
where M, K are the skew-symmetric matrices given by 
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The modified multisymplectic PDE can be used to derive the modified LECL and LMCL. 
Substituting z, M, K and S, into © and ©, the modified LECL and LMCL are found to be, 
respectively, 



E t + F x 
G x + I t 



At 2 
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0. 



where E, F, G, and / are given by equations (f9|)- (|10|) . In the next section, we numerically verify 
that these modified local conservation laws are satisfied to higher order. 



5 Numerical Results 

For our numerical experiments we consider the NLS equation with periodic boundary conditions, 
u(x + L,t) = u(x,t). We use initial data for a multi-phase quasi-periodic (in time) solution, i.e., 
uq(x) = 0.5(1 + 0.1 cos //a;), n = 2tt/L, L = 2\/2vr. This initial data corresponds to a multi-phase 
solution, near the plane wave, characterized by one excited mode. We examine the performance of 
the centered cell box scheme (which we designate as MS-CC) for varying mesh sizes and time steps. 
The solution to equation 1)24(1 is found by writing it in matrix form A^u 1 = A+u° + F{u l , u°) and 
using an iteration technique to solve for u . 

The solution with N = 64 and dt = 5 x 10~ 3 for 450 < t < 500 is shown in Figure The 
surface clearly exhibits the correct quasiperiodic behavior in time. In addition, we are interested 
in how well the local and global conservation laws are satisfied. To evaluate the local conservation 
laws, we use midpoint discretizations of the form 

1/2 _ fVg E 1/2 F 1 -F > 
V2 At + Ax ' 
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Figure 1: The MS scheme with N = 64 and t = 5 x 10 3 , T = 500: a) surface, b-c) error in the 
LECL and LMCL d-e) error in the global energy and momentum. 

In general, these residuals are not zero (Figures ^-c) . The errors in the local conservation laws, 
the LECL and LMCL ©-(jTUJ), are 10 -6 and 10 -3 , respectively, and are concentrated around the 
regions where there are steep gradients in the solution. If S(z) were a quadratic functional of z, 
S(z) = ^z T Az, with A a symmetric matrix, then the local conservation laws would be conserved 
exactly [2]. In general, as in the present case, the PDE is nonlinear and S(z) is not a quadratic 
functional. Therefore, the local energy and momentum conservation laws will not be preserved 
exactly. However the numerical experiments show that the local conservation laws are preserved 
very well over long times. In addition to resolving the LECL and LMCL very well, the MS scheme 
preserves the global errors extremely well. The error in the global energy oscillates in a bounded 
fashion, as expected of a symplectic integrator (Figure ^i) while the error in the global momentum 
(Figure^) and the norm (not shown) are conserved exactly (up to the error criterion of 10 -14 in 
the solver) since they are quadratic invariants. 

The maximum error in the LECL and LMCL and in the global energy and momentum for the 
MS scheme are provided in Table 1 for varying mesh sizes and time steps. The error in the LECL 
depends only on the timestep t and is second order, while the error in the LMCL depends only on 
the spatial mesh size N and is second order. 
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N 


32 


32 


32 


64 


64 


64 


t 


2.0E-02 


1.0E-02 


5.0E-03 


2.0E-02 


1.0E-02 


5.0E-03 


LE 


6.0E-05 


1.5E-05 


4.0E-06 


8.0E-05 


2.0E-05 


5.0E-06 


LM 


1.7E-02 


1.7E-02 


1.7E-02 


4.8E-03 


4.8E-03 


4.8E-03 


GE 


7.3E-05 


2.0E-05 


5.0E-06 


7.6E-05 


2.2E-05 


5.0E-06 


GM 


1.2E-13 


2.5E-14 


2.0E-13 


1.3E-13 


1.0E-13 


4.5E-13 



Table 1: The absolute maximum error in the local energy and momentum and the global energy 
and momentum obtained using the MS scheme MS, with T = 500. 

We next examine whether the modified local conservation laws obtained using the MS-CC 
discretization of the NLS are preserved to a higher order than the original local conservation laws. 
Since our solution is quasiperiodic, we compute the solution for < t < T, where T is chosen to 
include a characteristic cycle. From Figure^ T = 10 is sufficient. Since the ECL is independent 
of Ax (see Table 1), for a fixed N, we let At — ► 0. That is, start with Ato = Ax and let 
At = At /2 k , k = 0,1,..., 6. 

We compute the LECL and the modified LECL at each time step using centered approximations 
of the derivatives of sufficiently high order so as not to affect the order of the MS-CC discretization 
of the residuals. Figure © shows the loglog plot of the maximum error as a function of the timestep 
for the LECL and the modified LECL. Clearly we can see that while the LECL is satisfied to 2nd 
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Figure 2: Loglog plot of the error against At for the original ECL (o - -), and the MECL (x - -). 

order, the modified LECL is satisfied to 4th order. Verification of higher preservation of the LMCM 
as a function of the mesh size is similar. 
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